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We investigate theoretically how the stress propagation characteristics of granular materials 
evolve as they are subjected to increasing pressures, comparing the results of a two-dimensional 
scalar lattice model to those of a molecular dynamics simulation of slightly polydisperse discs. 

We characterize the statistical properties of the forces using the force histogram and a two-point 
spatial correlation function of the forces. For the lattice model, in the granular limit the force 
histogram has an exponential tail at large forces, while in the elastic regime the force histogram 
is much narrower and has a form that depends on the realization of disorder in the model. The 
behavior of the force histogram in the molecular dynamics simulations as the pressure is increased is 
very similar to that displayed by the lattice model. In contrast, the spatial correlations evolve qual- 
itatively differently in the lattice model and in the molecular dynamics simulations. For the lattice 
model, in the granular limit there are no in-plane stress-stress correlations, whereas in the molecular 
dynamics simulation significant in-plane correlations persist to the lowest pressures studied. 



I. INTRODUCTION 



Stress transmission in dry granular media is unusual because in these materials no simple relation exists between 
stress and strain |^-|^. Physical ingredients that give rise to this are that there are no tensile forces, that the 
particle deformations are very small, and that the particles can rearrange ||^ . Over the last several years evidence has 
accumulated that force propagation in dry granular media could be fundamentally different than in elastic solids 
||,|^-|^. Equations that have been proposed to describe stresses in lightly loaded granular media have the property 
that specification of boundary conditions at the top surface of the system is sufficient to determine the stresses 
throughout P,p|,|7|-pd|, in marked contrast to the elliptic equations of elasticity theory. 

However, applying a large enough uniform pressure to a granular material will cause it to exhibit an elastic linear 
response to a small additional stress. This is because uniform pressure both inhibits rearrangements (because it 
suppresses Reynolds dilatancy) and compresses the contacts, so that the non-tensile constraint on the interparticle 
forces becomes irrelevant. Thus, if stress propagation in lightly loaded granular media is indeed substantially different 
than in elastic media, then subjecting the material to high pressures will change fundamentally the stress propagation 
characteristics. 

This paper investigates theoretically the stress propagation in granular materials as they are subjected to increasing 
pressures. The goals of this work are to understand the physical mechanisms governing the evolution between granular 
and elastic behavior and to make specific experimental predictions for the behavior of granular media under increasing 
loads. 

We study a two-dimensional model system and compare the results to molecular dynamics (MD) simulations of 
two-dimensional systems of slightly polydisperse discs. Numerical studies of statistical models of granular media, 
where geometrical complexity is modeled in terms of uncorrelated random variables, are much faster and simpler than 
molecular dynamics simulations. Models of this type hold promise as a means to obtaining insight into the physics 
underlying the force propagation in granular materials. Our model for the granular regime is the two-dimensional scalar 
g- model Though the q- model has deficiencies ||l^, it is attractive because of its simplicity and its prediction of 

an exponential tail in the probability distribution of stress within a packing agrees with experiments p"3|-p^ and with 
simulations |l6| - |2l| . Our model for the elastic regime is a network of springs with a regular topology, with disorder 
introduced via randomly chosen spring constants p^-|2^]. To model the crossover between the two regimes, we exploit 
our observation that the q-model can be written as a scalar elastic network subject to certain constraints. Enforcing 
these constraints to an increasing degree, which causes the force propagation behavior to evolve from that of an elastic 
system to that of the g-model, models the crossover between elastic and granular behavior by a particulate assemblage 
subjected to decreasing pressure. 

We test the lattice model by comparing the results from the model to those of our MD simulations of two-dimensional 
systems of slightly polydisperse discs, focusing primarily on the probability distribution of stresses and on the two-point 
stress-stress correlation functions. The results of this investigation are mixed. The crossover in the force histogram 
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between the elastic network and the q-model is strikingly similar to the crossover observed in the molecular dynamics 
simulations as the pressure on the system is decreased. However, the lattice model and the molecular dynamics 
simulation exhibit qualitatively different trends in the behavior of the two-point correlation functions of the stress. 

The paper is organized as follows: Section |l| defines the scalar networks that we investigate. Section [D details 
the process of generation, solution, and analysis of these networks and discusses the generation of the molecular 
dynamics simulations of slightly polydisperse discs. Section |^ reports the results of the force distributions and 
spatial correlation functions for both the scalar lattice model and the MD simulations. Section ^ compares the results 
of the scalar lattice model, the MD simulations, and relevant experiments. Section VI summarizes and interprets our 
results. Appendix ^ calculates a finite-size correction to the in-plane stress-stress spatial correlation function for the 
g-model that is relevant to the interpretation of our numerical results. 



II. SCALAR ELASTIC NETWORKS AND THE Q-MODEL 



This section discusses the relationship between the g-model and the elastic network studied in this paper. Both 
models are scalar and are defined on a two-dimensional lattice. A scalar model is appropriate for a spring network 
if either the the network is very highly stretched |2^, ^ ^ , or if the motions are constrained so that displacements 
are unidirectional [ p5| . We consider the second situation and denote the direction along which the motion occurs as 
y, with positive y pointing downwards. 

Consider a network of nodes connected by springs on a diamond lattice as shown in Fig. ^, where the motion of 
every node is constrained to be along the vertical direction y. Each spring has the same unstretched length, so that 
in the limit of zero load the system forms a regular lattice. The springs connecting the nodes have spring constants 
that are chosen independently from a fixed probability distribution. Periodic boundary conditions are imposed in the 
horizontal direction, and the locations of the nodes at the top and bottom boundaries are fixed so that the vertical 
displacement of all the nodes in these rows relative to the unloaded configuration are identical. We index the nodes 
so that a node in column j in a row i with odd (even) i lies along the same vertical line as the other nodes in column 
j in rows with odd (even) indices. 

Let yij be the position of the node in row i and column j measured relative to its location in the absence of a load, 
and let fc' j and j be the spring constants of the springs emanating downward from the node at row i and column 
j. Every spring obeys Hooke's law, so that flj and f[j, the forces exerted on node (z, j) by the left and right springs 
below the node, are k\ .j{yi^j - yi+i,j-i) and klj{yi,j - Ui+ij) for odd i ikl j{yij - yi+i,j) and kljiytj - y^+i^j+i) for 
even i). The system is then compressed by setting yi j = [Ly — l)Ay and yLy,j = for all j, where rows 1 and Ly 
are the top and bottom rows, respectively, and AF is the average strain. We define Fi^j to be the total vertical force 
incident from above on node («, j), so that Fi^j = /JLij_i + fi-ij- The forces and displacements are determined by 
balancing the forces at every node, Fij = fl j + fij, and requiring that each yij be well-defined. This latter condition 
can be written as S = — [d]Y; here, S is the strain and Y is the displacement field [ p7[ . 

In our spring networks, each spring constant has a value selected independently at random from various probability 
distributions that are described below. We obtain the forces and strains along each link of each network using the 
method outlined in Ref. p7| ]. 

This scalar elastic model is equivalent to a resistor network [^|2^, Forces and strains in the elastic system 
correspond to currents and voltages, respectively, in the resistor network. The requirement that the vertical forces at 
each node balance is equivalent to Kirchhoff's current law, while the requirement that the position of each node is 
well-defined is equivalent to Kirchhoff's voltage law. 



A. Comparison between the elastic model and the g-model. 

The force Fij incident from above on node is transmitted to the sites below in the two pieces //^- and f[j. 
Because of force balance, one can always write 

4, =9^^.., /Z, = (1-9.,,)^;.,. (1) 

In a g-model, the g^.j are random variables that are chosen independently at every site. In an elastic network, Eq. (|l|) 
still holds, but the qij are determined by the configuration of random spring constants together with the requirement 
that the displacement field be single- valued. For spring constants that are chosen independently, the force along 
any branch will depend on the values of the spring constants throughout the system. Important consequences of 
this non-locality include the presence of spatial correlations between the (?ij 's and a nontrivial relation between the 
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distribution of spring constants and the distribution of the g's, including possibly the presence of g's that are negative, 
indicating the appearance of tensile forces in the network. 

A key observation underlying our work is that the g-model is equivalent to an elastic network subject to the 
constraint that the strain on every spring in each row is identical. The strain need not be constant from one row to 
the next, but it is simplest to consider the case in which it is. Let the amount of strain be Ay. Given the total force 
incident on node (i,j) from above, Fij, if one chooses the spring constants to be 



'^hj ~ ' ~ AY ' 



then the force exerted down the left link emanating from node is kl jAY — qijFij, and the force exerted down 
the right link from node i,j is k^jAY = (1 — qij)Fij. This force redistribution rule is exactly that of of the g-model. 
Given the set of qij values and the forces at each node in the top row of the system, we can create an equivalent 
spring network in a layer-by-layer manner. 

We do not implement explicitly a no-tensile force constraint in our networks, in contrast to the work of Refs. |2^] 
and However, in the g-model limit, there are no tensile forces. Our molecular dynamics simulations of lightly 
loaded material yield force distributions much closer to that of the g-model than to those of the non-tensile elastic 
networks of Ref . |^ . 

To study the crossover between elastic and g-model behavior, we generate iteratively a sequence of networks that 
interpolate between the elastic and q-model limits. The procedure adjusts the spring constants to make the strain in 
the system more uniform while keeping the ratio of spring constants emanating from each node constant. At iteration 
n, the spring constants klj{n) and k^ ^{n) are set to 

klAn)= ^-^^^^^\ l-q,,), (3b) 

where Fij(n — 1) is the force through node at iteration n — 1. The qij are kept fixed, and the iteration procedure 
is started with Fij{0) = 1. 

To characterize the crossover between elastic and q-model behavior as the iteration proceeds, we need to quantify 
the degree to which the constant-strain constraint is violated. We use as our measure of the spatial variation in the 
strain the dimensionless quantity 

Ly-i L 



6Sj, ^ (4) 

SY 

where SY^j = Yi^j - Kj+ij-i and 6Y^.j = Yij - Y^+ij for odd i {SY^^ = Y^j - Yi+ij and SY^.^ ^ Yi^j - K^+ij+i for 
even i), and 

^=MZ^^E>^^.+^^:.) = Ay (5) 

Here, Ly and are the number of rows and columns, respectively. In the elastic limit 8Sjq ~ 0.2, and as discussed 
above, 8Si^ is zero for the q-model. 



III. METHODS 



A. Scalar lattice model 

We consider diamond-shaped lattices with springs on each link, as shown in Fig. 0. The positions of the top and 
bottom node layers are fixed and periodic boundary conditions are imposed in the transverse direction. The forces 
along all the links depend on the choice of spring constants, {^ij}, and are calculated using the node-potential method 
described in Ref. The overall strain for each network is scaled so that the average vertical force through each 

node is normalized to unity. 
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(6) 



.1 j=i 



where the sum is over the nodes in the network. 

Networks of height Ly = 500 are used, with analysis performed on separate groups of layers to distinguish between 
edge and bulk effects. The widths ~ 16 and 128 are powers of 2 in order to take advantage of FFT techniques in 
the calculation of spatial correlation function values described below. The number of realizations averaged over varies 
from 10 to 50, depending on lattice size and number of iterations. 

For the elastic regime, we use four different distributions of spring constants: uniform distribution of k^^ for 
k^^ e (0, 1), gaussian distribution of k^^ with the configuration average {k~^) — 1 and standard deviation cTfe-i = 0.5, 
uniform distribution of k with k € (0,1), and gaussian distribution of k with fc = 1 and ak = 0.5. We construct 
networks with = 16 and 128 with 50 and 25 realizations, respectively. 

For the g-model regime, a uniform distribution of q with q e (0, 1) is used. We implement the iterative scheme with 
networks of size = 16 and 128 with 50 and 10 realizations, respectively, for 100 iterations. 

The local stress redistribution in a real granular material depends on microscopic details such as particle shape, 
friction characteristics, and preparation history. Instead of attempting to model the local force redistribution rules 
microscopically, our statistical models treat them as random variables chosen from different probability distributions. 
Since these probability distributions are not known a priori, we wish to identify and study properties that are not 
sensitive to the choice of the probability distribution governing the local force redistribution in the model. We focus 
on P{F), the probability distribution of stresses at the nodes; P{q), the probability distribution of the redistribution 
fractions q; and the spatial correlation functions of the force fluctuations about the mean values p8[. 
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where SFij^ — Fij — F and 6qij = qij — q; F is the average force and q is the average q value. The indices I and 
m in Eq. M label layers and columns, respectively, while k and j are the spatial separation in layers and columns. 
These correlation functions are normalized so that Co(0) = 1 and Co(0) = 1. Positive values (correlation) indicate a 
tendency for nodes separated by k rows vertically and j columns horizontally to be either both above or both below 
the mean, while negative values (anti-correlation) indicates opposite behavior of one above and one below the mean. 



B. Molecular dynamics simulations 



Here we discuss our molecular-dynamics (MD) simulations to used to generate 2-D packings of discs. Varying the 
ratio of external load pressure to particle stiffness induces crossover between granular and elastic behavior. We cal- 
culate the probability distributions and corresponding spatial correlation function values for forces and redistribution 
fractions q that are analogous to those in the scalar model. 

Our simulations employ a method similar to that used by Durian and collaborators for sheared foams [^9|-^, 
in addition incorporating kinetic friction, contact damping, and particle rotation, and using two different repulsive 
interparticle force laws (linear and Hertzian). 



1. MD interaction rules 



The discs in our simulation are all of identical mass mjj — 1 and interact via purely repulsive normal contact forces 
and kinetic friction. The interaction force between two discs whose centers are at positions f and fj with radii a.; and 
aj is non-zero only if their separation Sri_j < 0, where 

Sri^j = |fi - fjl - {ai + aj). (8) 

The normal contact force J^ij is calculated from the overlap I'Jrij I. We examine two force laws. The first is a linear 
force law based on a spring-like restoring force that yields 
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^i.j — Kij \Srij I , 



(9) 



with Kij — {1/Ki + l/Kj) ^, where Kd is the spring constant for disc d. The second is a non-hnear force law based 
on Hertzian contacts between spheres ^2| , 

+ (10) 
'-' at aj 

where D = | ((1 — cr^)/E), with a and E being the material properties Poisson's ratio and Young's modulus, re- 
spectively. For both force laws, the forces are directed so as to separate the overlapping discs. To calculate forces 
generated by interactions with walls, we assume the walls to be discs of infinite radius. 

Kinetic friction is incorporated into the disc interactions although static friction is not. The introduction of frictional 
forces causes the discs to rotate; however, the frictional force is zero at mechanical equilibrium. The kinetic friction 
force fij for contact between discs i and j is 

fzj = /i^jj, (11) 

where fi is the coefficient of kinetic friction, and is directed opposite to the contact point velocity if^^y For disc 
i, this velocity v^^j is related to disc velocities Vi and Vj, the angular velocities uji and ujj, and directional vector 
r ^ in- rj)/\ri - rj \ by 

= ~ ■'r)f + (flj'^i + ajLOj) X f, (12) 

where Vij — Vi — Vj. 

Damping during contact between discs i and j is used as an additional means of dissipating kinetic energy. It is 
generated by applying to disc i a force J-u and torque Fd given by 

= -Atrans^i, (13a) 
Td = -AangWi, (13b) 

where v'^ is the translational velocity of disc i relative to the interaction center of mass for the two discs i and j 
that are in contact and lo'^ is its angular velocity relative to the total angular momentum of the disc pair. Atrans and 
Aang are damping constants. This process conserves both translational and angular momentum. Energy is directly 
removed from the system as opposed to being converted between translational and rotational motion. 

The bottom and top walls have mass mw and are constrained to move only vertically. An inward force of magnitude 
-Fwaii is applied to each wall in order to compress the system. Damping of the wall motion suppresses volume oscillations 
and serves as the primary means of energy dissipation. The damping force J^wd on a wall is 

^WD — ~^wvw, (14) 
where vw is the velocity of the wall and X\y is the wall damping constant. 



2. MD implementation 

Ensembles of systems of iV = 1024 discs of average radius ajj are generated by starting with triangular array of ^/N 
rows and ^/N discs per row placed in a horizontally periodic system with both height and width L — 2.2TiaD\fN ■ For 
the data shown here, discs are placed in the system at positions {L{nx + 0.05) / a/ZV, L{ny + 0.5) / Vn) for odd Uy and 
{Lijix + 0.55) / \/N , L{ny + 0.5) / ViV) for even Uy, with indices Ux and Uy running from to ^/N — 1. In practice, discs 
with gaussian distributed polydispersity of a a ~ Q.lao placed on this triangular array do not overlap. The results 
obtained are not sensitive to initial disc placement. The system is then compressed by the application of an inward 
force on the top and bottom walls. All discs have the same spring constant Kd = K = 1. The coefficient describing 
wall damping is set to Xw/ttt-w — 1- Damping coefficients for translational and angular motion for disc contacts are 
set to Atrans/'Ti£) — 1 and Xangl Id — 4.1, where Id = ^''nDCijj is the moment of inertia for a disc with radius an- The 
coefficient of kinetic friction is set to /i = 0.2 for both disc-disc and disc-wall contacts. Comparisons with samples 
produced without disc-contact damping or kinetic friction revealed no measurable differences in force probability 
distributions or in the two-point force correlation function. Incorporating additional energy-dissipation mechanisms 
allows systems to reach mechanical equilibrium more rapidly. The end time for each compression stage is chosen so 
that the average residual kinetic energy for each disc is equivalent to translational movement of approximately or less 
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than O.OlaD in unit time. Because of the increased external energy input, systems at higher compressions are allowed 
a less restrictive limit of approximately O.OSau. Visual inspection of final configurations do not reveal significant 
fluctuations in time in contact network topology or force magnitude in load-bearing structures. Comparisons with 
test systems with longer run times also do not show any significant quantitative differences. 

For a system of fixed size L, the typical compression of the system can be controlled through variations in the disc 
spring constant K or applied external force i^waii- Typical relative particle deformations 5TZ is given by 



^-^11 (^) 



where Nc is the total number of contacts, the sums are over pairs of discs i and j in contact, and 11 = F^aii/i is the 
external pressure. This estimate is approximate due to geometric factors and distributional fluctuations; however, the 
scaling of deformations toIl/K should hold generally. In our simulations, the disc spring constant K is held fixed and 
the pressure 11 is varied to induce crossover between granular and elastic behavior. We define the reference pressure 
n = IIo such that the relative particle deformation 5TZ « 6.25 x 10""'. The reference compression pressure Ho yields 



a force histogram typical of the granular range, as discussed below in Sec. [V. After the initial compression with 
n = Hq, the applied pressure is increased in stages to 11 = lOOIIo, at which 6TZ ~ 0.01. We also decrease the pressure 
from the initial 11 = Ho configuration down to 11 = O.OlIIo {671 « 10"^) in order to approach the zero-deformation 
limit. Fig. shows a sample MD system subjected to the pressures O.lIIo, Ho, lOIIo, and 50no. 

2 

For spheres with Hertzian contacts (using Eq. [lO| ), the deformation can be approximated by i57?,[^^l « (^-^^^ ^ ■ 

For our simulations D is chosen to yield deformations of the same order of magnitude as the linear contacts at the 
compression 11 = Ho. The pressures studied are the same as for the linear spring contact systems. 



IV. RESULTS 



A. Scalar lattice model 



Here we present our results for the scalar lattice models. We study how the probability distribution of total vertical 
force F incident on a node from above P{F) and the two-point force correlation function Ck{j) characterize the 
behavior in scalar elastic lattice networks in which the constant-strain constraint is enforced to varying degrees. In 
the g-model both P{F) and in-plane force-force correlation function Co(j) exhibit robust behavior for generic choices 
of probability distributions of q's. We investigate the degree to which these quantities depend on the choice of spring 
constant distributions in the elastic networks, and discuss the crossover of P{F) and Co(i) between the elastic and 
q-model behavior as the constant-strain constraint is implemented with increasing accuracy. 



1. Results for the q-model 

In the q-model, the force histogram P{F) decays exponentially at large forces [0 and Co{j) is zero for non-zero 
j These properties hold for a wide variety of choices of the distribution of g- values. 

Our results for the crossover from elastic to q-model behavior are obtained for the specific choice that the g's are 
uniformly distributed in [0, 1]. A two-dimensional q-model with this distribution of g's yields 

P{F) = 4Fe-^^. (16) 

For a system of infinite lateral extent, the in-plane force-force correlation function Co(j) = Sjo where 6jo is the 
Kronecker delta function |ll[]. For a system of finite width L^, force correlations must arise because all forces are 
positive and the total force through a layer is fixed. As discussed in Appendix assuming that this mechanism is 
the only one giving rise to correlations, one obtains that a 2-D system of lateral extent has Co(j) given by 

C,{j ^ 0) = -(i. - (17) 

This form for Co(j) agrees with our numerical results for the q-model on lattices of finite width. 
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S. Elastic networks 



For elastic networks with different distributions of spring constants, the probabiUty distribution of vertical force 
P{F), shown in Figs. ^, is narrower than that of the g-model. Its functional form depends on the choice of spring 
constant distribution. Choosing the spring constants k from a distribution either uniform in k or gaussian in k yields 
P(F)'s that are roughly gaussian while the P(F)'s for networks for distributions uniformly distributed or gaussian in 
k~^ display a tail at large F that is consistent with an exponential decay. Networks with gaussian distributed k or 
k~^ exhibit narrower P(P)'s than their counterparts with uniformly distributed k or k~^. 

In contrast to the behavior of the force probability distribution P{F), the force- force correlation function values 
Cfc(j) are quantitatively indistinguishable for all the distributions of spring constants that we examined, as shown in 
Fig. ^. For Co(j), the force- force correlation function within the same layer, we see a strong anti-correlation for j = 1 
of magnitude ~ 0.30 that decays within j < 8. For vertical separation fc > 0, we see a simultaneous reduction in peak 
magnitude (at j = 0) and broadening of peak width but with the anti-correlation signature and decay joining the 
curve laid out by fc = 0. 

The probability distributions of redistribution fraction g, P{q), shown in Fig. are roughly gaussian and peaked 
a.t q = q = 0.5 for all distributions of spring constants examined. The widths of the P{q) depend on the choice 
of distribution of spring constants, with the gaussian distributed fc and k~^ once again being narrower (standard 
deviation aq « 0.16 and 0.15, respectively) than their uniformly distributed counterparts {aq ^ 0.25 for random k and 
(Tq ~ 0.21 for random fc~^). All of the elastic networks display significant correlations between g's at different nodes as 
demonstrated in Fig. ||, which shows the correlation function Ck{j) for all the random distributions. The correlations 
between g's are an important factor in determining the statistical distribution of the forces in these systems; Fig. |^ 
shows that a g-model system with the same P{q) as an elastic network with a uniform distribution of fc but with no 
correlations between the q's yields a P{F) markedly different from the elastic network. 

No differences between the bulk (layers 201-300) and edge (layers 1-100 and 401-500) sections are detected in the 
distributions P{F) and P{q) or the correlation functions Ck{j) and Ck{j)- The results for lattices with = 128 are 
the same within statistical errors to the results from = 16 lattices. 

In the elastic networks, forces in less than 1% of the branches are tensile, and no node in any of the networks is 
subject to a tensile net force. Our results for P{F) for uniformly distributed fc's are very similar to those reported in 
Sexton et al. ||25[| , where a non-tcnsilc force constraint is enforced. 

3. Iterated networks- the q-model limit 

We now discuss the networks generated by our iterative algorithm for converting an elastic network to a g-model 
system. First, we verify that the generated networks do eventually converge to the g-model. After 100 iterations, 
the forces along the links of the iterated spring network are identical to those of the corresponding g-model to within 

io-4_ 

A subtle point in the method is that our iterative scheme yields a configuration in which the forces at the top and 
bottom boundaries of the iterated network may have nonzero spatial correlations, as the initial iteration n = system 
is elastic. As one proceeds away from the top and bottom boundaries, these correlations decay via a diffusive process 
that takes on the order of layers ||ll|. Thus, forces at different sites in the same layer are effectively uncorrelated 
only for systems with large aspect ratios. This result is consistent with our numerical observation that in fully iterated 
systems, correlations between forces at different sites in the same layer are present throughout the = 128 systems, 
while they are only present in the top- and bottom-most 200 layers of ~ 16 systems. 

4- Iterated networks- crossover between elastic and q-model behavior 

In the iterated networks the target values of the qij are fixed at the outset of iteration procedure. The initial 
state (zeroth iteration) is an elastic network with spring constants given by fc- ^ = g^j/AF and klj = (1 — qij)/AY. 
Therefore, the initial probability of node forces P{F) and the spatial correlation function Co(j) are those of elastic 
networks with spring constants chosen from a uniform distribution of fc. The realized q distribution P{q) (as opposed 
to the distribution of the target q values) is peaked at g = 0.5 and its spatial correlation function Ck{j) reveals slight 
nearest- neighbor correlations for fc = and anti-correlations at j = for fc > 0, once again matching elastic-regime 
behavior. 

The probability distribution of node forces, P{F), is shown in Fig. ^(a) for different values of iteration number 
n. As the number of iterations is increased, the P{F) develops an exponential tail at large forces. Fig. ||(b) shows 
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the probability distribution of the g's, P{q)^ versus the number of iterations. P{q) approaches the target form of a 
uniform distribution after roughly 10 iterations. 

Fig. ^a) shows our results for nearest-neighbor in-plane and vertical force- force correlation function values Co(l) 
and C2 (0) as the number of iterations n is increased. Fig. |^(b) shows the corresponding force- fraction q-q correlations 
Cq{1) and 6*2(0). While only about 10 iterations are necessary before the nearest-neighbor spatial correlations between 
q values go to zero, in-plane force-force correlations are still present after 100 iterations although much reduced in 
magnitude from the initial elastic (iteration n = 0) value and approaching asymptotically the expected zero-correlation 
value. 

The quantity 5Sn that we use to characterize the crossover between elastic and g-model behavior is defined in 
Sec. ||. As Fig. |l^ demonstrates, we observe 5Sn decreasing as the number of iterations n is increased according to 
a power law. A fit that assumes the dependence is of the form 5Siq cx yields a = —1.68 ± 0.02. 



B. Results of molecular dynamics simulations 

Here we discuss the results of our MD simulations. Fig.|rT| shows P{F), the probability distribution of vertical forces 
F = T-y, ioT MD systems under various applied pressures H. As with the scalar model, F has been normalized so that 
the average vertical force F = 1 for each system configuration. The progression of P{F) as pressure is increased is very 
similar both qualitatively and quantitatively to the crossover from granular to elastic behavior in the scalar model 
lattice systems. We calculate the force- force correlation values Co(j), shown in Fig. ^(b), by defining discs to be in 
plane with a tolerance of iO.lOao and j in units of average disc diameter 2aD- In contrast with the scalar model 
behavior, the MD systems exhibit a significant nearest-neighbor anti-correlation for all applied external pressures. 
These results for P{F) and Co(j) are independent of whether the samples are compressed in stages or directly at a 
fixed pressure 11 . 

We define the q value of a disc as the fraction of total vertical force received from its topward neighbors that is 
transferred to its bottom leftward neighbors. The probability distribution of q values, P{q) is shown in Fig.|l^. We also 
calculate the q-q correlation values Co(j) and Ck{0), although the large errors prevent the extraction of quantitative 
trends. Narrowing the statistical errors would be computationally prohibitive. 

The number of contacts increases significantly with the pressure, as shown in Fig. ^ As the magnitude of the 
typical overlap increases, additional contacts are formed. The number of contacts at low pressures is below the 
theoretically predicted average of Z = 2(i |^ , where d is the dimension of the system, because the polydispersity in 
radii and the lack of gravity allow for the existence of "rattlers" which do not support any of the external load. 

Our results for the Hertzian contact systems are indistinguishable from those of the linear springs throughout most 
of the range of pressures explored. At higher pressures H (H > SSHq), the added stiffness of the Hertzian contacts 
leads to the slower narrowing of P{F). 



V. COMPARISON OF RESULTS OF MD SIMULATIONS AND OF SCALAR ELASTIC NETWORKS 

Here we compare the behavior observed in the MD simulations and in the scalar elastic networks. Because different 
schemes are used to induce the granular-elastic crossover in the two systems (iterations in the scalar networks and 
external pressure for MD), we need to establish a common measure to quantify a system's position within the crossover 
region. As the evolution of the probability distribution of vertical forces, P{F), is qualitatively and quantitatively 
similar in the network model and in the MD simulations, we use matches in its form to establish a relationship between 
iteration number n and applied pressure H. Fig. |lj(a) shows matches in form between linear-force-law MD packings 
and iterated scalar network systems for H/Hq = 100 and iteration n — 0, H/Hq = 10 and n = 10, and H/Hq = 1 
and n = 100 systems. From these matchings, we map the iteration number n in the scalar networks to the equivalent 
applied pressure Heq(n) in the MD using the simple scaling 

noq(») ^ ioo_ ^^g^ 

Ho n 

We perform a check on this proposed scaling by considering the analogous quantities of deviation from constant 
strain in scalar systems SSn, given by Eq. ^ and deviation from the infinitely hard, zero-deformation limit in MD 
systems calculated by 



Nc f-^, [a^ + ajY 
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where Nc is the total number of contacts and the sum is over pairs of discs i and j in contact. We match 5S values for 
n/IIo = 10 and n = 10 by scaling the square deviation for the scalar network systems by a constant factor of 0.030. 
Fig. [lj(b) shows that this scaling yields reasonable agreement between 5Sn and 5Smd over the crossover region. 
In contrast to the agreement in the trends of P{F), qualitative differences exist between the scalar network model 



and the MD simulations in spatial correlation function values Cj{k). Fig. 15 [a) shows the nearest-neighbor in-plane 
and vertical force-force correlation values, Cq{1) and C2(0), for the crossover between elastic and granular regimes. 
While the MD systems exhibit a significant in-plane nearest-neighbor anti-correlation throughout the crossover, a 
decrease in its magnitude is seen in the scalar networks as the systems change from elastic to granular. MD systems 
do not exhibit strong vertical correlations, in contrast with the scalar networks whose C2(0) value increases significantly 
as the granular limit is approached. 

The large statistical uncertainties in our q-q correlation functions for MD systems restrict us to making only 
qualitative behavior descriptions. The trend for in-plane nearest-neighbor correlation behavior Co(l) in both systems 
is similar. However, qualitative differences exist for vertical correlation value C2(0): the MD systems display consistent 
anti-correlation behavior, while the scalar networks display anti-correlation behavior in the elastic regime which decays 
rapidly to uncorrelated behavior as the granular limit is approached. 

Our work indicates that experiments on granular media at high pressures should yield a force histogram that differs 
qualitatively from that observed at lower pressures. Experiments by Howell et al. [ p6| as well as experiments and 
simulations by Makse et al. |^ are in qualitative agreement with this result. Howell et al. control the transition 
between granular and elastic behavior of slowly sheared systems in a 2-D Couette geometry by varying the packing 
fraction 7 within a range 0.77 < 7 < 0.81. The average force/length on a particle increases with 7. For lower values of 
7, the distribution of large stresses is asymptotically exponential, while the distribution of stresses has a gaussian form 
at higher packing fractions 7. Makse et al. pTf apply increasing pressure to three-dimensional packings of spherical 
glass beads to achieve the crossover between granular and elastic behavior and also perform MD simulations on 3-D 
systems. Makse et al. observe a crossover in the force histogram P{F) in a pressure range that is consistent with our 
2-D MD results. 

An interesting question is whether the persistent in-plane nearest-neighbor anti-correlation in the forces that is ob- 
served in the MD simulations is present in experimental systems. Mueth et al. [ [Tst do not find evidence of correlations 
between different sites in the same horizontal layer; any nearest-neighbor anti-correlation in the experiment is smaller 
than the experimental resolution. However, they measure a different correlation function, Ki{r), defined as 

Nb Nb 

E E ^('^^^ - 

i—l 

where the sums are over the Nb particles in the bottom layer, fi is the force at position in the bottom layer and 
rij ^ \ri — rj\. Calculation of Ki{r) from the numerical data for our MD simulations yields values of the correlation 
function that are smaller than the error bars in the experiment. Comparison with Ref. is necessarily qualitative 
since the experiments measure the properties at the surface of a 3-D packing while our MD results are calculated 
using numerical data from the bulk of a 2-D system. 



VI. DISCUSSION 



We have investigated the crossover between elastic and granular stress transmission in both a 2-D scalar lattice 
model and in molecular dynamics simulations of slightly polydisperse discs. The evolution of P{F), the probability 
distribution of stresses, is very similar in the lattice model and in the MD. However, the behavior of the spatial 
correlation functions for stress, Cfc(j), differs qualitatively. 

Our investigations of the scalar model have several implications for the development of granular media models. First, 
we have shown that implementing a local constraint can convert an elastic network to a q-model. This constraint has 
the natural physical interpretation that the strain in the system must be uniform; it is plausible that rearrangements 
would prevent strain gradients from forming. Second, implementing this constraint to increasing accuracy causes 
the force histogram P{F) to evolve in a manner similar to that observed in the MD simulation as the pressure is 
decreased. P{F) has a tail consistent with exponential decay at large forces in the granular limit, while the P{F) for 
the highly compressed system is much narrower and decays more quickly at large forces. We note that implementing 
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a non-tensile force constraint alone, as in Ref. |25|, yields gaussian decay in P{F) at large forces even at the lowest 
pressures, in qualitative disagreement with the MD results of us and others [p^|l8|,^. 

While this success in describing the evolution of the force histogram and the scalar model's simplicity in both 
formulation and implementation make it an attractive platform for the study of media models, the discrepancy in the 
behavior of the correlation function behavior with the MD simulation results needs to be addressed. The scalar model 
assumes explicitly that in the granular regime the stress redistribution fractions q at different sites are uncorrelated. 
The extent to which this condition is valid needs to be examined in more detail. Spatial correlations of the g's can 
strongly affect the probability distribution of stress P{F) but the degree to which these correlations exist in 

real packings has not been settled. A possible source of spatial correlations in the g's is the constraint that non-tensile 
vector forces must be balanced. However, vector generalization of g-model systems that have been proposed to date 
have required arbitrary constraints to be imposed to limit the scale of stress components perpendicular to the direction 
of applied force ||l2| . Clarification of the roles of vector force balance and contact formation is key to identifying and 
characterizing the processes governing stress transmission beyond those that have been implemented in the scalar 
model. 

In conclusion, we have shown that similarities exist in the evolution of the probability distribution of stresses 
P{F) in the crossover between elastic and granular regimes for a scalar lattice model and MD simulations of slightly 
polydisperse discs. However, the systems exhibit qualitative differences in the two-point force correlation function 
Ck{j). Further investigation of the systematic influences leading to the spatial correlations between forces is necessary 
for the development of a successful model of stress transmission in granular media. 
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Sid Nagel, Joshua Socolar, and Bob Behringer for useful conversations. This research was supported by the MRSEC 
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APPENDIX A: FINITE-SIZE CORRECTION TO CORRELATION FUNCTION CALCULATION 

In the g-model in the limit of infinite size, forces at different sites in the same layer are completely uncorrelated. In 
a system of finite transverse extent, the requirements that the total force through every layer is identical and there are 
no tensile forces lead to a finite-size correction to the correlation function. This appendix discusses this correction. 

We characterize the correlations between force fluctuations on sites in the same row using the correlation function 



where SFi^m — Fi,m — F is the deviation of the force at a site in row I and column m from the average force F. For 
the g-model with a uniform distribution of g's, in a system of infinite transverse extent this correlation function is |Tl| ] 



This result follows from the fact that P(Fq,,F^), the probability that the force through node a is Fa and the force 
through node /3 at the same horizontal level is F^, is factorizable: 
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On a lattice of finite width [L^ sites) , thc^ multipoint force probability distribution function must be consistent with 
the facts that first, the total force down every layer is fixed, and second, no force is negative. This implies that 

• The maximum force on any node in any layer cannot be larger than Fmax = L^F, and 

• The force Fa at a node a contributes to the total force along a layer and hence affects the sum of the forces 
through the remaining sites in the layer. Defining F as the average force through all the sites in the layer other 
than site a, we have 
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Assuming that the only correlations present in the finite system are those required to satisfy these conditions, the 
joint probability distribution in the system with finite can again be written 



nF.,F,) = j^p(^)p{^), 



but now the distributions are subject to the constraints 



P{uj)duj =1, / LuP{uj)duj = 1 

10 Jo 

/ P{u)di^ = 1, / 

Jo Jo 



vP{v)dv = 1 . 



(A5) 

(A6a) 
(A6b) 



In the limit of oo, we expect corrections to P{F) to be of order 1/Lx as the fluctuations in the forces at the 

sites arc of order unity. As we will see, with the assumptions that we have made, the finite size correction to the 
correlation function does not depend on the form of the probability distribution P. Taking the new constraints into 

account, for a ^ P we have 
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As the correlation function is normalized with respect to the average fluctuation size, 
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which is independent of the form of the P{F). As oo, Co{j 7^ 0) ^ 0, as expected. 
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FIG. 1. Elastic network considered in this paper. Each node is connected to two neighbors in the row above and to two 
neighbors in the row below, with movement confined to the vertical direction. A node and its surrounding nodes have 

been labeled. The system is compressed by holding the top and bottom rows each at fixed positions. Disorder in the stress 
distribution is introduced by variation of individual spring constants k. In the elastic regime, the k values are chosen at random 
while in the granular regime, the fc's are additionally constrained so that the strain in each row is constant. 
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(a) n = o.ino 



(b) n = Ho 




(c) n = lOHo (d) n = 5ono 

FIG. 2. Contact networks of a sample MD-genoratod packing of 1024 discs for different values of applied external pressure 
n. The reference pressure Hq is such that the average fractional change in particle radius at a contact is 6.25 x 10"'*. In the 
transition from granular to elastic behavior, the number of contacts in the system increases and the magnitudes of the contact 
forces become more homogeneous. While the width of contacts shown is a proportional to the force magnitude, they have been 
rescaled for each pressure so direct comparisons between subfigures is not possible. 
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FIG. 3. Force probability distribution P{F) for 128-column random spring configurations. The form of P{F) in the elastic 
regime depends on the distribution of spring constant values k. However, these P{F) for all distributions of k occupy a narrow 
range in comparison with the q- model granular regime result shown by the solid gray line, (a)-(d) We fit the functional form 
P{F) oc F^e~^^ to the random k~^ distributions with A = 6.67 and B = 7.95 for the uniform distribution and A = 14.70 
and B — 16.21 for the gaussian distribution. For the random k distributions, we fit P{F) oc e'^~^' with 5" = 0.47 for the 
uniform distribution and S = 0.32 for the gaussian distribution. No differences are seen between layer groups near the edge or 
in the bulk of the systems. 
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FIG. 4. Spatial force-force correlation function Ck{j) for the forces within the bulk layer grouping (layers 201-300) in 
128-column elastic- regime scalar networks. In-plane correlations Co{j) are shown in the main plot while the vertical correlation 
peak decay given by Cfe(O) is shown in the inset. An observed nearest-neighbor in-plane anti-correlation appears to be robust 
with respect to variations in spring constant distributions. The vertical correlation is similarly robust. 
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FIG. 5. Force fraction probability distribution within the bulk layer grouping (layers 201-300) for 128-column elas- 

tic-regime scalar networks. The fits are gaussians peaked at q = 0.5 with the width being dependent on the spring constant 
distributions. Gaussian distributed and fc configurations are narrower, with widths a ~ 0.16 and 0.15, respectively, than 
their uniformly distributed counterparts (random k, a ~ 0.25 and random fc~^, a ~ 0.21). 
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FIG. 6. The two-point spatial correlation function Ck{j) for the force fraction q within the bulk layer grouping (layers 
201-300) in 128-column elastic-regime scalar networks. In-plane correlations Co{j) are shown in the main plot with vertical 
correlations Cfc(O) in the inset. As with the force-force correlations, varying the spring constant distribution has minimal effect 
on these correlations. 
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FIG. 7. Effect of spatial correlation of q values on the probability distribution of force P{F). Force fraction probability 
distributions P{q) of the uniform distribution of k elastic network and a g-model system generated by choosing q values from a 
gaussian distribution centered at g = 0.5 with width cr, = 0.25 is shown in the inset. While the two P{q) distributions appear 
nearly identical, the spatial correlations in the elastic network yield a functionally different form for the probability distribution 
of forces P{F) than that of the uncorrelated g-model system. 
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FIG. 8. (a) Force probability distribution P{F) at various stages of iteration for the bottom layer grouping (401-500) of an 
ensemble of 16-column scalar networks. Initially P{F) is similar to the distribution for uniformly distributed k systems, shown 
by the grey dashed line. The distribution broadens with increasing iterations with small F agreement with granular g-model 
systems being achieved on order of 20 iterations. Further iterations are necessary to approach agreement for large values of F. 
The P{F) distribution for g-model systems is shown by the black dashed line, (b) Force fraction probability distribution P(q) 
at corresponding stages of iteration. The distribution of q values approaches the expected uniform distribution within the first 
10 iterations. 
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FIG. 9. (a) Nearest-neighbor in-plane and vertical force-force correlation values Ca(l) and C2(0) at various stages of iteration 
for an ensemble of 16-column scalar networks. The observed anti-correlation in nearest-neighbor forces decreases in magnitude 
as the number of iterations increases and appears to approach asymptotically the expected zero-correlation value. The vertical 
correlation increases in magnitude with increasing iterations, indicating a stronger preference for the formation of force vertical 
channels, (b) Nearest-neighbor in-plane and vertical q-q correlation values Co{l) and (72(0). The spatial correlations for the 
force fraction q decrease rapidly in magnitude as their values approach the uncorrelated target q distribution. 
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FIG. 10. Deviation from constant strain 5Sn for each layer in a scalar network system defined by Eq. 0. The solid line is 
a power law fit of form 5Sn oc n~°' with a = —1.68 ± 0.02. Increasing iterations confirm the approach to the constant strain 
limit, which is equivalent to the q model. 
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FIG. 11. (a) Progression of vertical force probability distribution P{F) between elastic and granular regimes in MD-generated 
packings. At high applied pressure 11, P{F) has a gaussian functional form, as shown by the gray dashed line, similar to the 
elastic behavior of the scalar network. At low pressure, P{F) has widened and displays an exponential tail. For comparison, 
P{F) for the g-model is shown by the black dashed line. The apparent transition to the granular regime appears to occur by 
roughly II/llo — 1. (b) In-plane force correlation function value Co(j) for various applied pressures. We define discs to be 
in-plane within a tolerance of iO.lOao with spacings j given in units of average disc diameter 2aD- Unlike the scalar networks, 
in the MD we see a consistent nearest- neighbor ant i- correlation for all pressures. 
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FIG. 12. Progression of the vertical force fraction q distribution P{q) between elastic and granular regimes in MD-generated 
packings. Each q is calculated as the fraction of the vertical force component that is transferred to the bottom left neighbors 
of a disc. At high pressures H, we see a peaked form centered at g = 0.5. For decreasing pressure, P{q) flattens out and is 
roughly uniform. The increasing magnitude of the leftward bin with decreasing pressure If is due to increasing probability of 
isolated "rattlers" and discs with no bottom left neighbors. 
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FIG. 13. Number of contacts Nc and coordination number Z for MD-generated packings of 1024 discs versus externally 
applied pressure 11. The number of contacts fits roughly to a form Nc ~ No + a(n/no)'', where A^o is the number of contacts 
in the zero-force hmit. At low pressures, A'o is slightly less than two contacts per particle because of the presence of "rattlers" 
in the zero-gravity system. 
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FIG. 14. (a) Matching vertical force probability distributions P{F) for iterated scalar spring networks and MD packings of 
particles with a linear force law. We observe good agreement in the form of P{F) between iteration n — and II/Ilo — 100, 
n = 10 and Il/llo = 10, n = 100 and ll/llo — 1. (b) Normalized strain deviation SSn and SSmd as a function of actual 
or equivalent applied pressure. The matching of P{F) is used to generate a mapping between iteration values of the scalar 
networks and the pressures imposed on the MD systems. A simple scaling approach yields llcq/llo — 100/n. Additionally, SSn 
for the scalar networks is scaled by a factor of 0.030 to achieve equivalence at n = 10 (n/IIo — 10). 
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(a) Force correlations - MD and scalar networks 
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FIG. 15. (a) Comparison of nearest-neighbor in-plane and vertical force-force correlation function values Co(l) and C2(0) 
between MD and scalar spring networks at various pressures 11 and iterations n. We see qualitative and quantitative differences 
in behavior. MD systems exhibit a significant in-plane nearest-neighbor anti-correlation throughout the crossover region. This 
contrasts with the decrease seen in the scalar networks as it approaches the granular limit (large iteration n). Strong vertical 
correlations develop in the scalar networks but not in the MD systems, (b) Nearest-neighbor in-plane and vertical q-q correlation 
function values Co{l) and (72(0). The sizeable errors in the MD values allows for only qualitative comparisons. Both systems 
exhibit similar trends for in-plane nearest-neighbor correlation values. However, qualitative differences exist for the behavior 
of vertical correlations: the MD systems display a constant anti-correlation throughout the crossover region while the scalar 
networks exhibit an anti-correlation in the elastic regime which decays rapidly to uncorrelated behavior as the granular limit 
is approached. 
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